part 2

## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpuQYs4g/downloaded_packages

## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpuQYs4g/downloaded_packages

## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpuQYs4g/downloaded_packages

load(here("jk_code", "JK_cleanMD.rds"))

1 Removing mt, Rp, Gm genes

SO4 <- subset(SO4, features = grep("^mt-|^Rp|^Gm", rownames(SO4), invert = TRUE, value = TRUE))

ElbowPlot(SO4)

SO4 <- SCTransform(SO4) %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:15) %>%
    FindClusters(resolution = 0.20) %>%
    RunUMAP(dims = 1:15)

## Running SCTransform on assay: RNA

## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.

## Calculating cell attributes from input UMI matrix: log_umi

## Variance stabilizing transformation of count matrix of size 15822 by 11431

## Model formula is y ~ log_umi

## Get Negative Binomial regression parameters per gene

## Using 2000 genes, 5000 cells

## There are 1 estimated thetas smaller than 1e-07 - will be set to 1e-07

## Found 286 outliers - those will be ignored in fitting/regularization step

## Second step: Get residuals using fitted parameters for 15822 genes

## Computing corrected count matrix for 15822 genes

## Calculating gene attributes

## Wall clock passed: Time difference of 16.30294 secs

## Determine variable features

## Centering data matrix

## Place corrected count matrix in counts slot

## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Set default assay to SCT

## PC_ 1 
## Positive:  Umod, Egf, Tmem52b, Sult1d1, Sostdc1, Fabp3, Klk1, Prdx5, Foxq1, Cox6c 
##     Wfdc15b, Ly6a, Wfdc2, Krt7, Slc25a5, Ckb, Cox5b, Atp5g1, Cldn19, Cox4i1 
##     Ggt1, Atp1a1, Ndufa4, Cox8a, Atp5b, Chchd10, Gadd45g, Cox6b1, Atp1b1, Cox7b 
## Negative:  Pappa2, Zfand5, Aard, Robo2, Pde10a, Wwc2, Itga4, Nadk2, Nos1, S100g 
##     Ramp3, Neat1, Sgms2, Ptgs2, Col4a4, Irx1, Col4a3, Tmem158, Mir6236, Ranbp3l 
##     Itprid2, Bmp3, Cdkn1c, Camk2d, Sdc4, Mcub, Dctd, Etnk1, Srrm2, Peg3 
## PC_ 2 
## Positive:  Fos, Junb, Jun, Hspa1b, Hspa1a, Btg2, Egr1, Zfp36, Atf3, Ier2 
##     Fosb, Socs3, Klf6, Klf2, Jund, Dnajb1, Dusp1, Tsc22d1, Gadd45g, Rhob 
##     H3f3b, Gadd45b, Cebpd, Ubc, Mt1, Mt2, Actb, H2bc4, H1f2, Klf4 
## Negative:  Mir6236, Pappa2, CT010467.1, Egf, Slc12a1, Umod, Etnk1, Wnk1, Sfrp1, Atp1b1 
##     Nme7, Robo2, Sptbn1, Col4a4, Hsp90b1, Sec14l1, Pde10a, Zbtb20, Itprid2, Dst 
##     Mal, App, Lars2, Wwc2, Col4a3, Kcnq1ot1, Atrx, Utrn, Tfcp2l1, Pou3f3 
## PC_ 3 
## Positive:  Fth1, Ubb, Ftl1, Ldhb, Car15, Fxyd2, Prdx1, Cd63, Cd9, Eif1 
##     Mpc2, Mgst1, Mt1, Clu, Bsg, Tmem213, Aldoa, Mdh1, Ramp3, Itm2b 
##     Spp1, S100a1, Selenow, Tmem59, Tmem176b, Wfdc2, Tmbim6, Tspo, Atpif1, Ppp1r1a 
## Negative:  Mir6236, Egf, CT010467.1, Umod, Neat1, Tmem52b, Nme7, Fos, Malat1, Slc12a1 
##     Kcnq1ot1, Jun, Junb, Lars2, Etnk1, Dst, Egr1, Wnk1, Slc5a3, Sult1d1 
##     Foxq1, Atrx, Atp1b1, Pnisr, Fosb, Atf3, Ivns1abp, Zfp36, Btg2, Syne2 
## PC_ 4 
## Positive:  Pappa2, Aard, Tmem52b, Umod, Egf, Sult1d1, Tmem158, Foxq1, Mcub, Ptgs2 
##     Dctd, Wwc2, Iyd, Car15, Ramp3, Ptprz1, Hsp90b1, Cd9, S100g, Wnk1 
##     Cdkn1c, Defb42, 1700028P14Rik, Pth1r, 5330417C22Rik, Tmsb4x, Tagln2, Ctsc, Clu, Bmp2 
## Negative:  Mt1, Apoe, Mt2, Aebp1, Sostdc1, Gpx6, Fxyd2, Ptger3, Neat1, Ckb 
##     Fgf9, Egfl6, Ivns1abp, Mfsd4a, Defb1, Car4, Plet1, CT010467.1, Fkbp11, Chchd10 
##     Atp5md, 2900052N01Rik, Cox6c, Igfbp5, Atp5k, Tmem213, Chka, Mgst3, Avpr1a, Atp1a1 
## PC_ 5 
## Positive:  Hspa1a, Hspa1b, CT010467.1, Pappa2, Klk1, Car15, Cldn10, Fth1, Hspa8, Jun 
##     Lars2, Fau, Mir6236, Itm2b, Wfdc2, Uba52-ps, Aard, Ftl1, Ptger3, Eef1a1 
##     Pik3r1, Egf, Sfrp1, Mal, Gpx4, Tspo, Atp1a1, Tmem176a, Id3, Hsp90aa1 
## Negative:  S100g, Actb, Tmem52b, Sdc4, Tmsb10, Ppia, Abhd2, Uroc1, Egr1, Serf2 
##     Slc39a1, Igfbp7, Atf3, Ndufb1-ps, Atp5md, Alkbh5, Cebpd, Cox6c, Ramp3, Ndufa3 
##     Gnb1, Atp5e, Foxq1, Atp5k, Ldhb-ps, Fam107a, Kdm6b, Rbm47, Atp5mpl, Wfdc15b

## Computing nearest neighbor graph

## Computing SNN

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 11431
## Number of edges: 355780
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8800
## Number of communities: 5
## Elapsed time: 1 seconds

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 12:06:11 UMAP embedding parameters a = 0.9922 b = 1.112

## 12:06:11 Read 11431 rows and found 15 numeric columns

## 12:06:11 Using Annoy for neighbor search, n_neighbors = 30

## 12:06:11 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 12:06:11 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpuQYs4g/file1181a258d4435
## 12:06:11 Searching Annoy index using 1 thread, search_k = 3000
## 12:06:13 Annoy recall = 100%
## 12:06:13 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 12:06:14 Initializing from normalized Laplacian + noise (using RSpectra)
## 12:06:14 Commencing optimization for 200 epochs, with 480882 positive edges
## 12:06:14 Using rng type: pcg
## 12:06:16 Optimization finished

DimPlot(SO4)

DimPlot(SO4,split.by = "sample")

DimPlot(SO4,split.by = "treatment")

SO4@meta.data <- SO4@meta.data %>% 
  mutate(subclass_MD = dplyr::case_when(
    seurat_clusters == 0  ~ "type_1",
    seurat_clusters == 1  ~ "type_1",
    seurat_clusters == 2  ~ "type_2",
    seurat_clusters == 3  ~ "type_3",
    seurat_clusters == 4  ~ "type_4",
  

  ))

SO4@meta.data$subclass_MD <- factor(SO4@meta.data$subclass_MD , levels = c("type_1", "type_2", "type_3", "type_4"))

Idents(SO4) <- SO4@meta.data$subclass_MD

DimPlot(object = SO4, reduction = "umap", group.by = "subclass_MD", label = TRUE)

DimPlot(object = SO4, reduction = "umap", label = TRUE)

Idents(SO4) <- "subclass_MD"

DimPlot(SO4,split.by = "sample",group.by = "seurat_clusters")

SO.markers <- FindAllMarkers(SO4, only.pos = TRUE)

## Calculating cluster type_1

## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session

## Calculating cluster type_2

## Calculating cluster type_3

## Calculating cluster type_4

SO.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)

SO.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10
DoHeatmap(SO4, features = top10$gene) + NoLegend()

## Warning in DoHeatmap(SO4, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay: Ifi47,
## Rsad2, Gbp9, Gbp3, Oasl2

type_1_markers <-(SO.markers[SO.markers$cluster == "type_1", ])
type_2_markers <-(SO.markers[SO.markers$cluster == "type_2", ])
type_3_markers <-(SO.markers[SO.markers$cluster == "type_3", ])
type_4_markers <-(SO.markers[SO.markers$cluster == "type_4", ])

#Type 1

df<- type_1_markers %>% arrange(desc(avg_log2FC))

df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% rownames_to_column(var="SYMBOL")

markers <- DEG_list %>% 
  rownames_to_column(var = "SYMBOL") 


ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)

## 'select()' returned 1:1 mapping between keys and columns

## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 1.79% of input gene IDs are fail to map...

markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC))) %>% 
arrange(p_val_adj) %>%
head(100)
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go

## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:100] "22682" "56089" "268902" "68646" "23850" "23984" "52357" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...288 enriched terms found
## 'data.frame':    288 obs. of  12 variables:
##  $ ID            : chr  "GO:0001822" "GO:0072001" "GO:0072073" "GO:0051926" ...
##  $ Description   : chr  "kidney development" "renal system development" "kidney epithelium development" "negative regulation of calcium ion transport" ...
##  $ GeneRatio     : chr  "13/97" "13/97" "8/97" "6/97" ...
##  $ BgRatio       : chr  "375/28928" "390/28928" "173/28928" "76/28928" ...
##  $ RichFactor    : num  0.0347 0.0333 0.0462 0.0789 0.0519 ...
##  $ FoldEnrichment: num  10.34 9.94 13.79 23.54 15.46 ...
##  $ zScore        : num  10.56 10.31 9.79 11.41 9.77 ...
##  $ pvalue        : num  4.15e-10 6.68e-10 1.28e-07 2.20e-07 3.73e-07 ...
##  $ p.adjust      : num  7.80e-07 7.80e-07 9.94e-05 1.28e-04 1.60e-04 ...
##  $ qvalue        : num  5.54e-07 5.54e-07 7.06e-05 9.11e-05 1.14e-04 ...
##  $ geneID        : chr  "Robo2/Epcam/Irx1/Col4a4/Col4a3/Cdkn1c/Sdc4/Irx2/Enpp1/Bcl2/Hoxb7/Nrp1/Iqgap1" "Robo2/Epcam/Irx1/Col4a4/Col4a3/Cdkn1c/Sdc4/Irx2/Enpp1/Bcl2/Hoxb7/Nrp1/Iqgap1" "Robo2/Epcam/Irx1/Sdc4/Irx2/Bcl2/Hoxb7/Iqgap1" "Nos1/Mcub/Ptgs2/Bcl2/Calm1/Calm2" ...
##  $ Count         : int  13 13 8 6 7 9 7 9 6 6 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722

 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2 lines

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

2 Type 2

  df<- type_2_markers %>% arrange(desc(avg_log2FC))


df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% 
  rownames_to_column(var = "SYMBOL") 

ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)

## 'select()' returned 1:1 mapping between keys and columns

## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 2.94% of input gene IDs are fail to map...

markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC))) %>% 
  arrange(p_val_adj) %>% head(50)
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go

## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:50] "330428" "53315" "15220" "69824" "105243" "13645" "242653" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...81 enriched terms found
## 'data.frame':    81 obs. of  12 variables:
##  $ ID            : chr  "GO:0006119" "GO:0009060" "GO:0019646" "GO:0045333" ...
##  $ Description   : chr  "oxidative phosphorylation" "aerobic respiration" "aerobic electron transport chain" "cellular respiration" ...
##  $ GeneRatio     : chr  "19/50" "20/50" "15/50" "20/50" ...
##  $ BgRatio       : chr  "154/28928" "206/28928" "66/28928" "271/28928" ...
##  $ RichFactor    : num  0.1234 0.0971 0.2273 0.0738 0.1852 ...
##  $ FoldEnrichment: num  71.4 56.2 131.5 42.7 107.1 ...
##  $ zScore        : num  36.4 33.1 44.2 28.7 39.8 ...
##  $ pvalue        : num  5.26e-31 1.71e-30 9.00e-29 4.89e-28 2.68e-27 ...
##  $ p.adjust      : num  4.72e-28 7.68e-28 2.69e-26 1.10e-25 4.82e-25 ...
##  $ qvalue        : num  3.34e-28 5.43e-28 1.90e-26 7.75e-26 3.41e-25 ...
##  $ geneID        : chr  "Cox5b/Cox6c/Uqcrq/Chchd10/Uqcr10/Cox5a/Uqcr11/Ndufb8/Cox6b1/Cox7a2/Cox4i1/Cox7c/Ndufb9/Uqcrh/Cox8a/Cox6a1/Cox7b/Ndufa13/Uqcrb" "Idh2/Cox5b/Cox6c/Uqcrq/Chchd10/Uqcr10/Cox5a/Uqcr11/Ndufb8/Cox6b1/Cox7a2/Cox4i1/Cox7c/Ndufb9/Uqcrh/Cox8a/Cox6a1/"| __truncated__ "Cox5b/Uqcrq/Uqcr10/Cox5a/Uqcr11/Ndufb8/Cox7a2/Cox4i1/Cox7c/Ndufb9/Uqcrh/Cox8a/Cox6a1/Cox7b/Uqcrb" "Idh2/Cox5b/Cox6c/Uqcrq/Chchd10/Uqcr10/Cox5a/Uqcr11/Ndufb8/Cox6b1/Cox7a2/Cox4i1/Cox7c/Ndufb9/Uqcrh/Cox8a/Cox6a1/"| __truncated__ ...
##  $ Count         : int  19 20 15 20 15 15 20 15 15 8 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722

 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2 li

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

3 Type 3

  df<- type_3_markers %>% arrange(desc(avg_log2FC))


df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% 
  rownames_to_column(var = "SYMBOL")

ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)

## 'select()' returned 1:1 mapping between keys and columns

## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 24.26% of input gene IDs are fail to map...

markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC))) %>% 
  arrange(p_val_adj) %>% head(50)
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go

## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:50] "14282" "13653" "12702" "11910" "14281" "16598" "16477" "22695" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...446 enriched terms found
## 'data.frame':    446 obs. of  12 variables:
##  $ ID            : chr  "GO:1903706" "GO:0009408" "GO:0042026" "GO:0009266" ...
##  $ Description   : chr  "regulation of hemopoiesis" "response to heat" "protein refolding" "response to temperature stimulus" ...
##  $ GeneRatio     : chr  "12/48" "7/48" "5/48" "7/48" ...
##  $ BgRatio       : chr  "458/28928" "96/28928" "24/28928" "177/28928" ...
##  $ RichFactor    : num  0.0262 0.0729 0.2083 0.0395 0.0571 ...
##  $ FoldEnrichment: num  15.8 43.9 125.6 23.8 34.4 ...
##  $ zScore        : num  13 17.2 24.9 12.4 14 ...
##  $ pvalue        : num  8.94e-12 2.34e-10 4.21e-10 1.70e-08 2.15e-08 ...
##  $ p.adjust      : num  1.46e-08 1.91e-07 2.30e-07 6.96e-06 7.03e-06 ...
##  $ qvalue        : num  7.59e-09 9.92e-08 1.19e-07 3.60e-06 3.64e-06 ...
##  $ geneID        : chr  "Fos/Zfp36/Jun/Hspa1b/Tsc22d1/Zfp36l1/Hsp90aa1/Snai2/Nfkbiz/Irf1/Actb/Hspb1" "Hspa1a/Hspa1b/Dnajb1/Hsp90aa1/Ier5/Hspb1/Dnaja1" "Hspa1a/Hspa1b/Hsp90aa1/Hspb1/Dnaja1" "Hspa1a/Hspa1b/Dnajb1/Hsp90aa1/Ier5/Hspb1/Dnaja1" ...
##  $ Count         : int  12 7 5 7 6 9 8 7 7 8 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722

 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2 li

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

4 Type 4

  df<- type_4_markers %>% arrange(desc(avg_log2FC))


df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% 
  rownames_to_column(var = "SYMBOL")


ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)

## 'select()' returned 1:1 mapping between keys and columns

## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 39.47% of input gene IDs are fail to map...

markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC))) %>% 
  arrange(p_val_adj) %>% head(50)
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go

## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:23] "15945" "15957" "23962" "55932" "626578" "236573" "58185" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...59 enriched terms found
## 'data.frame':    59 obs. of  12 variables:
##  $ ID            : chr  "GO:0035456" "GO:0035458" "GO:0009615" "GO:0051607" ...
##  $ Description   : chr  "response to interferon-beta" "cellular response to interferon-beta" "response to virus" "defense response to virus" ...
##  $ GeneRatio     : chr  "11/22" "10/22" "13/22" "12/22" ...
##  $ BgRatio       : chr  "77/28928" "68/28928" "412/28928" "329/28928" ...
##  $ RichFactor    : num  0.1429 0.1471 0.0316 0.0365 0.186 ...
##  $ FoldEnrichment: num  187.8 193.4 41.5 48 244.6 ...
##  $ zScore        : num  45.3 43.8 22.8 23.6 44.1 ...
##  $ pvalue        : num  1.55e-23 1.63e-21 3.64e-19 2.24e-18 3.76e-18 ...
##  $ p.adjust      : num  5.66e-21 2.97e-19 4.43e-17 2.04e-16 2.74e-16 ...
##  $ qvalue        : num  3.23e-21 1.70e-19 2.53e-17 1.17e-16 1.57e-16 ...
##  $ geneID        : chr  "Ifit1/Gbp3/Tgtp2/Ifi47/Iigp1/Tgtp1/Ifit3/Gbp2/F830016B08Rik/Gbp7/Ifitm3" "Ifit1/Gbp3/Tgtp2/Ifi47/Iigp1/Tgtp1/Ifit3/Gbp2/F830016B08Rik/Gbp7" "Cxcl10/Ifit1/Oasl2/Rsad2/Tgtp1/Nlrc5/Oasl1/Rtp4/Ifit3b/Ifit3/Gbp2/Gbp7/Ifitm3" "Cxcl10/Ifit1/Oasl2/Rsad2/Nlrc5/Oasl1/Rtp4/Ifit3b/Ifit3/Gbp2/Gbp7/Ifitm3" ...
##  $ Count         : int  11 10 13 12 8 8 8 6 4 5 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722

 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2 li

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

save(SO4, file = here("jk_code", "JK_remove_mtrpgm_part2.rds"))

5 Make Heat map for DEGs that display function of each cell type

gene_list <- c("Ptgs2", "Nos1", "Mcub", "Egf",
              "Jag1", "Sulf2", "Plau",
              "Fosb", "Fos", "Junb", "Jun",
              "Cxcl10")

DoHeatmap(SO4, features = gene_list) + NoLegend()